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We consider the dynamics of the disordered, one-dimensional, symmetric zero range process in 
which a particle from an occupied site k hops to its nearest neighbour with a quenched rate w{k). 
These rates are chosen randomly from the probability distribution f{w) ~ (w — c)", where c is 
the lower cutoff. For n > 0, this model is known to exhibit a phase transition in the steady state 
from a low density phase with a finite number of particles at each site to a high density aggregate 
phase in which the site with the lowest hopping rate supports an infinite number of particles. In 
the latter case, it is interesting to ask how the system locates the site with globally minimum rate. 
We use an argument based on local equilibrium, supported by Monte Carlo simulations, to describe 
the approach to the steady state. We find that at large enough time, the mass transport in the 
regions with a smooth density profile is described by a diffusion equation with site-dependent rates, 
while the isolated points where the mass distribution is singular act as the boundaries of these 
regions. Our argument implies that the relaxation time scales with the system size L as with 
z — 2 + l/{n + 1) for n > 1 and suggests a different behaviour for n < 1. 

I. INTRODUCTION 

The presence of quenched disorder is known to strongly affect the dynamical and steady state properties of many 
systems. For instance, for noninteracting particles moving in a random medium, disorder can lead to anomalous 
transport and change the manner in which the steady state is approached Moreover, when inter-particle 

interactions are present, disorder can induce new collective effects such as phase separation Q. In this paper, we 
study the effect of quenched, sitewise disorder on a simple stochastic model of interacting particles known as the zero 
range process [^D- This process can be viewed as describing a system of interacting particles hopping in and out of 
wells with various depths. While the static properties of this model are known analytically, the temporal properties 
are not as well characterized. Here we study the approach to the steady state of a zero range process which undergoes 
a disorder-induced phase transition. 

The study of stochastically evolving lattice models has played a central role in better understanding interacting, 
statistical systems [Q. Such models are defined directly through simple stochastic rules, in contrast to the traditional 
route of defining the interactions through a Hamiltonian and then constructing dynamical rules consistent with it. 
The zero range process describes one such class of lattice models; other important classes include the asymmetric 
simple exclusion process which is a simple model of a current-carrying system and the contact process whose 

steady state shows a phase transition from an active phase to a dead phase where no further evolution is possible [|| . 

The zero range process deals with a conserved number of particles hopping on a lattice with, in general, site- 
dependent rates. For almost all choices of these rates, the particles have on-site interactions; however, these rates 
do not depend on the particle occupation at other sites. In this sense, the range of interaction is zero. The steady 
state of this process is known to be of product measure form in all dimensions In the recent past, the zero range 
process has been used to model several systems with quenched disorder, such as traffic flow in a system of cars with 
different preferred speeds |9) and activated flow down a rugged slope besides several other applications A 
recent application is in a model of polymerization in a random medium with imperfect traps ||ll|. Interestingly, for 
several choices of disorder, as the density is increased, the steady state of the system shows a phase transition from 
a homogeneous phase in which the density is roughly uniform, to an infinite aggregate phase in which the density 
' profile has a singularity at the site with the lowest hopping rate. 

The latter state is particularly interesting since in this phase, starting from an initial random distribution of 
particles, an infinite aggregate is formed at the site with the lowest hopping rate. The question arises: what is the 
mechanism by which the system locates the site with globally minimum rate and transports a finite fraction of the 
total number of particles to it through the random medium? In this paper, we address this question for a symmetric 
zero range process in one dimension. Figure |^ shows how the density profile evolves in time starting from an initial 
random distribution of particles (Fig. ^) to the steady state in which an aggregate builds up at the site with the 
lowest hopping rate (Fig. Illd). 
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The rest of the paper is organized as follows. We define the model and discuss its steady state properties in 
Section ||. In Section III, we first give a simple, qualitative picture of the dynamics of relaxation to the steady state. 



The scaling properties of the relaxation time define a dynamic exponent which is determined using local equilibrium 
arguments, supported by Monte Carlo simulations. Finally, we conclude with a discussion of open questions. 

II. THE MODEL AND ITS STEADY STATE 

In this section, we define the model and briefly discuss its steady state properties. We consider the unbiased or 
biased motion of a conserved number M of particles on a d dimensional lattice of length L with periodic boundary 
conditions. At any site k occupied by a nonzero number of particles, a single particle attempts to hop out at a rate 
w{k), independent of the number of particles present at site k or its neighbours. Note that this choice of rates implies 
an attractive on-site interaction since the hop-out rate of these particles is lower than that of noninteracting particles. 
The zero range process with this particular choice of rates has appeared in several contexts, as in a particlewise 
disordered, asymmetric exclusion process [^|2| and as a limiting case of a model of aggregation and fragmentation 
with sitewise disorder [11 . 

As discussed in Section J, the problem is essentially to understand how the system locates the site with the minimum 
hopping rate. For this reason, it is useful to assign two labels to the hopping rate. We denote the rates by Wj{k) 
where k is the site index and j is the index when the rates are arranged in ascending order with j = 1 labelling the 
lowest rate. Correspondingly, rnj{k) denotes the mass (or number of particles) at the site with hopping rate Wj{k) 
at a particular instant. In the following, we will only display the labels which are pertinent to the discussion and 
suppress the rest. Also, unless the time dependence is explicitly specified, all the quantities will refer to the steady 
state. 

The rates {w} are chosen independently for all k = 1, L"^ from a common probability distribution 

f{w)^[{n + l)/{l-c)'"+^] (w-c)", we[c,l] , OO , n>0 . (1) 

For both the symmetric and the totally asymmetric case, the probability P{{'m{k)}) of a configuration 
{m{l),m{2), ...,m{L'^)} is known to be given by [D 

where Af is the normalisation constant and v is the fugacity which can be determined using the conservation law 

Sfc=i "^e^) = This solution holds in all dimensions and for all bias. In the presence of bias, Eq.(^) describes 
a nonequilibrium steady state whereas in the absence of bias, it describes an equilibrium state as the condition of 
detailed balance holds. This equilibrium state is described by a Hamiltonian which is long ranged, thus allowing the 
system to have a phase transition even in one dimension The probability P{m, k) that there are m particles at 
site k is given by 



which implies that the average number of particles at site k is 

w{k) — V 1 — s[k) 

where s{k) = J2m=i ^i^'^) ~ v/w{k) is the average occupation probability of site k. Since the total number of 
particles is conserved, 

P^J2 + / /H , (5) 

where p = M/L'^ is the total particle density and disorder averaging has been done. This equation is reminiscent 
of Bose-Einstein condensation in the ideal Bose gas where the lowest momentum state is macroscopically occupied 
beyond a critical density pc |12[. At p = pc, the fugacity v gets pinned to the lower cutoff of w, namely c (to which 
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wi will tend in the thermodynamic limit) so that the site with the lowest hopping rate supports an aggregate with a 
finite fraction of total particles. The condition p — 1(c) = determines the critical point pc = c(n + l)/n(l — c) where 
I{v) is the integral on the right hand side of Eq.(||). 

For p < Pc, the typical mass at all sites is of order unity. As the density is increased, there is a phase transition 
a,t p — Pc- In the high density phase, the slowest site supports mass of 0{L'^) whereas the site with rate Wj, j ^ 1 
supports mass of 0((L'^/j)i/"+i). The latter can be seen via a simple argument. Consider a variable x distributed 
uniformly between and 1. Since on average, L'^ observations of x will be equally spaced, it follows that the jth 
lowest observation in L'^ trials is typically at a distance j / L"^ above zero. Since (toj) is inversely proportional to the 
separation between the lowest and the jth lowest hopping rate and f{w) can be related to the uniform distribution 
by a change of variables (w = c + (1 — c) x^^^~^^), one obtains {rrij) ^ [L"^ / j)^^"^^^ , j ^ 1. 



III. APPROACH TO THE STEADY STATE 



In this section, we will discuss the approach to the steady state for the symmetric, zero range process in one 
dimension with the hopping rates chosen from j(w) defined in Eq.(^. Wc first describe a simple picture of approach 
to the steady state. Then we present an argument based on local equilibrium which suggests that the relaxation time 
scales with system size L as U' with z = 2 + l/(n + 1) for n > 1. However, this treatment breaks down for n <\. 

We first illustrate qualitatively the temporal sequence of events through which an aggregate with mass of 0(V) is 
formed at the site with the lowest hopping rate, starting from a random initial condition in which each site has mass 
of order unity (Fig. |l|a). The behaviour of the average mass (mj(t)) for j = 1, ...,4 as a function of time is shown 
in Fig. ||. Here (...) is to be understood as an ensemble average over evolution histories. We find that {mi{t)) rises 
steadily and then saturates to its steady state value while each of {m2{t)), {1713(1)) and {m^lt)) rise, attain a maximum 
and then decay to their respective steady state values. This non-monotonic behavior is not hard to understand. At 
some finite time, each particle is able to move only a finite distance away from its initial position and tends to get 
temporarily trapped at the site with the lowest w within the neighbourhood explored by it, rather than the global 
minimum wi . A typical configuration at such intermediate time thus has several large aggregates at such local minima 
while the rest of the system has masses of order unity (Fig. |l|b). This explains why at short enough times, {mj(t)) for 
j = 1, ...,4 are of the same order and increasing (Fig. ^). As time progresses, the particles are able to access larger 
regions in space and identify new local minima. Then the mass increases at these newly accessed local minima at the 
expense of the previous ones (Fig. ^). This explains the drop in {mi{t)) and (m3(t)) after they have reached their 
respective peak values, though (mi(t)} and {m2{t)) continue to rise (Fig. |^). Finally, {m2{t)) also starts dropping 
and the excess mass is transported to the location of the global minimum (Fig. |^d). Once the global minimum is 
recognised by all the particles, the system reaches a steady state and {mj{t)) (wj) for all j. 

We now turn to an analytical description of the mechanism of the mass transport. The exact time evolution equation 
obeyed by {m{k,t)) for a given {w} can be written as 

^^^^^^ = w(/fc-l)s(fc-l,i) + w(A: + l)s(fc + l,i)-2w(fc)s(fc,i) . (6) 

At large enough times, the system is expected to be in local equilibrium. This allows one to assume the s{k, t) — {m{k, t)) 
relation to be approximately as in the steady state (Eq.(U)). Substituting for s(fc, t) in the above equation, one obtains 

^^^^^^^=G(fc-l,i) + G(fe + l,t)-2G(fc,i) , (7) 

where G(fc, t) = w{k){m{k, t)) /(I + (m(fc, t))). 

The treatment so far is valid for any density, but we are primarily interested in densities for which the system is 
in the infinite aggregate phase in the steady state. For such densities, at large times, one needs to divide the system 
into two sets - set A composed of those sites which support rather large aggregates and set B at which the masses 
are small and close to their steady state values (see Fig. ||). The elements of these sets are not fixed in time; with the 
passage of time, the number of elements in set A reduce and that of set B increase by the mechanism described at the 
beginning of this section. Eventually, when the system is close to the steady state, set A is left with a few isolated 
sites that support aggregates whose masses scale as a nonzero power of L while the bulk of the sites belong to set B. 

We are interested in the fluctuations in the density profile about the steady state, Am(A;,t) = {m{k,t)) — {m{k)). 
Since Am is small for k B and very large for k G A, an expansion in Am{k,t) would be justified only for the 
background sites. Below we first analyse these background {B) sites and then treat the isolated sites in set A as 
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the boundaries of the B regions. For the B sites, we may retain the lowest order terms in an expansion of Eq.(^ in 
powers of Am and obtain 

dAm{k, t) ^ ^ _ ^ - 1, t) + D{k + 1) Am(fc + 1, t) - 2D{k) Am(k, t) , keB , (8) 
ot 

with 

D{k)^{w{k)~c)^/w{k) . (9) 

Here we have used the fact that v tends to c in the infinite aggregate phase. 

Equation (j^) describes a random walker in a random medium with site-dependent hopping rate D{k). Using Eqs.(|^) 
and we find that these rates are distributed according to the probability distribution g{D) ~ for D ^ 0. 

Note that whereas f{w) has a nonzero lower cutoff, g{D) has zero as the lower limit due to which it diverges as D ^ 
for n < 1. The problem of a particle moving with random, spatially inhomogeneous hopping rates is well studied 
and has been reviewed in For a configuration of randomly distributed {D} in a large system observed on large 
enough time scales, the mean squared displacement of the random walker grows as Vt where the diffusion constant 
V = 1/{{1/D{k))) and ((...)) stands for disorder average, provided {{1/D{k))) is finite Thus for large spatial 
and time separations, the particle diffusion is described by a single, effective diffusion constant. In our case, one can 
calculate {{1/D)) by noting that {{1/D)) = dl/dv\v^c- Expanding I{v) for v close to c, we obtain 

I{v) ^ Pc + 0{c-v) , n>l (10) 
= p, + 0{{c-vr) , 0<n<l , (11) 

which implies that {{1/D)) is finite for n > 1 and diverges for n < 1. The divergence in the latter case indicates 
anomalous diffusion i.e. the mean squared displacement grows sublinearly. In the remainder, we will restrict ourselves 



to n > 1 |13| 



On averaging over disorder configurations, for a system of size L, the mass {{mi{t, L))) on the site with lowest 
hopping rate is expected to scale as 

{{m,{t,L)))^tf'H(±) (12) 



where 



xj/ N [ constant for a; ^ 1 /-,o\ 
^^"^^Ix-'^ for.»l ■ (^^^ 



Since {{mi{t, L))) ^ L in the steady state, it follows that (3z = 1. We now determine the relaxation time of the 
system by estimating the time T required for {mi{t, L)) to reach its steady state value for a typical realisation 
of disorder. At large enough times, the background sites in set B form the bulk of the system whereas only a few 
isolated sites with low hopping rates belong to set A. Let us now consider time scales above which the two sites with 
the lowest rates (i.e. wi and W2) are the only elements left in set A. As shown in Fig. ^ (m2(i)) drops after it has 
reached its peak while {mi{t)) keeps rising at the expense of the former so that there is a net transfer of mass from 
the site with j = 2 to the site with j = 1- One can think of these two sites as the boundaries to the bulk background 
region with the former site feeding particles at rate W2 to the bulk and the latter at rate wi. On large time scales, a 
quasi-equilibrium is established, and in view of the discussion above, the particles diffuse through the bulk with an 
effective diffusion constant P. Thus there is a transfer of mass at a rate 'DAwi2/ri2 where Awi2 — W2 — wi and ri2 is 
the spatial separation between the two sites. The relaxation time T12, which is essentially the time taken to transfer 
the peak mass m2 at the site with rate W2, is given by 

_ ri2m2 

^-^''-VA^2 ■ ^ ^ 

Let us estimate the size dependence of the quantities that appear on the right hand side of the above equation. 
Typically the separation ri2 is of order L in which case m2 is also of order L. This is because for times earlier than 
when m2 peaks, the site j = 2 is the slowest one encountered by the particles in a finite fraction of the system. 
Consequently, the numerator in the above equation is proportional to L^. This result holds even in the exceptional 
case when the sites with j = 1 and j = 2 are separated by a distance of order unity, as they behave effectively as a 
single slow site and one can use the same reasoning as above on replacing the site with j = 2 by that with j = 3. 



4 



Further, using the argument given at the end of Section ||, we find that the inverse rate separation {Awi2)~^ scales 
as Collecting all the dependences, it follows that the relaxation time T scales with system size L asT 

where 

z = 2 + l/(n+l) , n>l . (15) 

We measured the growth of {{mi{t, L))) using Monte Carlo simulations and find that it grows as . As shown in 
Fig. 1^, our expression for P = 1/z — {n + l)/(2n + 3) is consistent with the numerics. 

We conclude this paper with a discussion of two open questions: (i) We have seen that {{1/D)) diverges for n < 1 
indicating anomalous diffusion, which suggests that the expression for z in Eq.(|l5|) may be invalid for n < 1. It would 
be interesting to find how the system relaxes in this case, (m) For the case with asymmetric hopping rates, using 
the arguments analogous to those given above, we would expect the dynamic exponent to be z = 1 + l/(n + 1) for 
n > 1. Once again, we may anticipate a different result for n < 1. In |^ and fl^ the dynamics of the asymmetric 
version has been analysed using numerical simulations and extremal statistics arguments for a related deterministic 
model. The expression for the dynamic exponent in 14| is the same as quoted above but their arguments do not 



make a distinction between the regimes below and above n — 1. The elucidation of the anomalous regime for both 
the symmetric and asymmetric zero range process remains an interesting open problem. 

We will not attempt a summary of this paper. Professor N. Kumar once pointed out that since it is a part of the 
whole, a proper summary should include a summary of the summary and a summary of that summary, and so on. 
We would rather not try ! We are very pleased that this paper will appear in the Festschrift volume for Prof. Kumar, 
whose own work has brought out many of the surprises that disordered systems have to offer. 
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FIG. 1. Density profile {m{k,t)) vs. k in the aggregate phase for a given realization of disorder at time 
(a) t = 2, (fe) t = 2^°,{c) t = 2^^,(d) t = 2^". The site with the lowest hopping rate is located at fc = 1 and the 
second lowest at fc = L/2. At t = 2^", the system is close to the steady state and (mi(l,t)) = 426. Parameters used: 
L = 256,p = 4,n = 2,c = 0.5. 
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FIG. 2. Time dependence of the average density {mj(t)) for four slowest sites namely j = 1 (squares), 2 (triangles), 3 (circles) 
and 4 (pentagons) for a given disorder configuration. The horizontal lines denote the steady state value for {rrijit)) calculated 
using Eqs.(|4|) and (H). Parameters used: L = 256, p = 4, n = 2, c = 0.5. 
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FIG. 3. log-log plot of {{m.i{t,L))) vs. t to show the growth exponent /3. The data for {{mi{t, L))) has been averaged over 
50 histories and 21 disorder configurations. The theoretical prediction for /? = 0.428 for n = 2 is plotted as a solid line for 
comparison. Parameters used: L = 16384, p = 4, c = 0.5. 
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